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SUMMARY 



A modification of Hansen f s lunar theory is given in a form that 
permits a purely numerical treatment of solar perturbations of 
planetary satellites or of perturbations caused by the earth in the 
motion of hypothetical lunar satellites. The development of the co- 
ordinates and of the velocities of a satellite is obtained in the form 
of a trigonometric series in four arguments with numerical coef- 
ficients. At the earlier stages of the computation we are not com- 
pelled to use Hansen's basic W-function: We make use of the fic- 

titious mean anomaly and the replacement "bar" operation — two 
characteristic features of Hansen's theory— to set up a process of 
iteration in a convenient form, rather than to compute w. The pref- 
erence given to the method of iteration was motivated by the cir- 
cumstance that input information must not necessarily be limited to 
the choice of the rotating ellipse as an intermediary orbit. A sys- 
tematic investigation of the orbital stability of lunar satellites can 
be based on the development given herein. 
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ON A MODIFICATION OF HANSEN'S 
LUNAR THEORY* 

by 

Peter Musen 

Goddard Space Flight Center 


INTRODUCTION 

In this paper a modification of Hansen's lunar theory (Reference 1) is given in a form that per- 
mits a purely numerical treatment of solar perturbations of planetary satellites as well as the per- 
turbations caused by the earth in the motion of hypothetical lunar satellites. The development of 
perturbations of very close lunar satellites does not present any difficulty and therefore is not given 
here. The case of more distant satellites was treated by Lemekhova (Reference 2) on the basis of 
Delaunay's theory for orbits with low inclinations toward the orbit of the moon. At the present time 
we are not in possession of a complete analytical theory of satellites with highly inclined orbits, al- 
though Delaunay's theory can be modified to become valid for such orbits. Because of this circum- 
stance, the use of a semi-analytical theory is suggested, as given in this paper. 

Hansen's original lunar theory requires the use of the derivatives of the disturbing functions 
with respect to the eccentricity and the parameters determining the orbit plane's position. This 
means that a literal development of the disturbing function must be obtained before the numerical 
values of the elements can be substituted. However, with the advent of electronic machines, more 
emphasis can be put on the development of purely numerical theories. 

Unlike Hansen’s lunar theory, his planetary theory (Reference 3) is a purely numerical one. For 
this reason the theory of perturbations in the orbit plane, as presented herein, takes Hansen’s plane- 
tary theory, and not his lunar theory, as a pattern. In addition, at the earlier stages of the computa- 
tions we are not compelled to use Hansen's basic w -function. Instead, we make use of the fictitious 
mean anomaly and the replacement TT bar" operation to set up the process of iteration in a convenient 
form, rather than to compute w. The preference given to the method of iteration over the method of 
computing perturbations of different orders was motivated by the circumstance that input information 
must not necessarily be limited to the choice of Hansen's rotating ellipse as an intermediary orbit. 

In addition, the programming and the computational process become more uniform if the method of 
iteration is used. 


•Also to be published in the Journal of Geophysical Research. 
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Only after the process of iteration is completed is the W-function formed in order to determine 
the perturbations in the mean anomaly and in the radius vector. 

The orbit of the disturbing body is taken to be a moving ellipse. 

The inclination of the satellite's orbit plane toward the orbit plane of the disturbing body can be 
large. For this reason we do not deduce the standard development of the latitude into series, but 
make use of four "redundant" parameters introduced by the author in the theory of artificial satellites 
(Reference 4). The introduction of these parameters makes all basic arguments linear in time from 
the outset, and it permits us to deduce the rectangular components of the position and velocity vec- 
tors very easily. 

The knowledge of these components facilitates the computation of osculating elements. 

To simplify the form of the development, we take — as did Hansen — the orbit plane of the dis- 
turbing body as a basic reference plane. 

The effect of this plane's motion on the relative motion of the satellite is very small and can 
easily be taken into consideration. 

The form of the disturbing function in our case favors only the choice of time as an independent 
variable, but not the true or the eccentric anomaly as in the case of an artificial satellite of the earth. 


DEVELOPMENT OF THE DISTURBING FUNCTION 

For a satellite of negligible mass, we have 



(See Appendix A for list of notations.) Putting 



(i) 


( 2 ) 


2 


we deduce: 


and 


a a 
■l - ' r 7 


3 2 1 2 

T s "TP 


ad 

r Tr 


an 


m' a 3 

t 

a 

a' 4 

7^ 

m ' a 4 

t 

a 

a' 5 

r 

m # a 2 

/ 

a 

a' 3 

7 

m ' 

X 'J 

af 

+ a 

a 

'4 

m ' 

4- /I 

a 4 


5 3 

2 s 


!p ! ») 


4 15 3 3 A 

-o-s - - 4 -s 2 p 2 + g-p 4 ) 


3 2 
2 s 


TP 2 ) 


5 3 1 2 

2 s 2 p s 


m'a 4 a' /35 15 , 3 \ 

4 a' 5 ‘ TT ‘ \ 8 s " T s P + gP j 


( 3 ) 


(4) 


a r ' dQ 
r (9 S 


= 3 


mla^ . a* /IS 2 3 2 

a' 4 ^ \ 2 S “ 2P 


•' a 4 a' /35 . 15 , \ 

a ' 5 ■ T 7- ^T s " T"P s j 


(5) 


We have 


s - cos ( v — v> ) cos + Sin (v“^) sin (v'-^/') cos I , 


r a 

a • yr cos ( v iyC/ ) co s { v * ~ 0 1 ) 


( 6 ) 


_r a 

+ a * 7 T sin (v - i/i) sin { v' ) cos I . 


( 6 ') 


Taking the relations 


2N - i// Q + ^ 0 ' “ P “ P* “ 2n 0 at 

2K = “ ’A/ ~ ^ ^ + 2 n 0 77 1 


3 


v “ f + tt q + n 0 yt , 
v' = f' + ^ + n 0 y' t 

into account, we deduce from the Equations 6 and 6' that 
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The ai'guments w and a/ are defined as linear functions of time from the very start. 
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These five series represent an input. They can be easily computed, even using a desk calculator, 
by means of Cayley's scheme (Reference 5). We have: 
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In this scheme we multiply every line of the upper part by every line of the lower part and add 
the arguments together, then attach the sin or cos symbol to each term of the product. The upper 
sign of the argument in the solar part is associated with the development of 
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PERTURBATIONS IN THE ORBIT PLANE 


To determine the perturbations in the radius vector and in the mean anomaly, Hansen from the 
very start introduces the function 
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From the standpoint of the process of iteration the use of the elements 
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would be preferable. We have (Reference 1) 
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we deduce from Equation 16: 
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Equation 14 takes a very simple form: 
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Equations 18, 21, and 24 bear more resemblance to the equations from Hansen’s planetary theory 
than to the corresponding equations from his lunar theory, and they can be used to serve as a basis 
of a purely numerical development. The value of n 0 y must be determined in such a way that no 
purely secular term is present in Equation 21. 


9 



PERTURBATIONS OF THE POSITION OF THE ORBIT PLANE 

In previous work (Reference 4) the author has established the following equations for the deter- 
mination of a parameters: 
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Substituting Equation 32 into Equations 27 to 30 and taking Equations 33 to 36 into account, we deduce 
the system of equations that replace Hansen's equations for his P, Q, and K parameters: 
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If necessary, the effect of the secular motion of the disturbing body's orbit can be easily taken into 
account. We have for this effect (Reference 6): 
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We can neglect the effect of change of i 7 on the motion of the satellite. Substituting 
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we obtain the effect of the secular motion of the disturbing body's orbit plane on the relative position 
of the satellite’s orbit plane in the form: 
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The combined effect of direct perturbations with the secular motion of the orbital plane of the dis- 
turbing body is given by the system of equations representing the sum of Equations 37 to 40 and of 
44 to 47. 
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The values of a and 77 must be determined in such a way that no terms linear in time are present in 
Equations 49 and 50. 
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DETERMINATION OF CONSTANTS OF INTEGRATION 


The series for ¥ has the form 
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and no additive constant of integration appears in this development. The secular motion n Q y of the 
perigee of the satellite in its orbit plane is determined in such a way that no constant term is present 
in Equation 21 for d'P/dt and, consequently, no secular term appears in Equation 52. The series for 
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or in the form 
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The value of n 0 iz taken from the previous approximation can be used in the computation of the right- 
hand sides of Equations 63 to 66 for our purposes. 
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The generalized Hill's formula (Reference 7) 
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The value of B is computed on the basis of the previous approximation. Designating by [w 0 J the value 
of W 0 as obtained using the formal integration procedure, we have 
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Two constants of integration, c, and c 2 , must be adjusted in such a way that no secular term and no 
term with the argument g is present in n 0 Sz. The part 
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contains the term of the form 


A t + A 2 cos g ■ 


Let P be the coefficient of the term in ( pf a 0 ) cos qp + ( 3/2) e 0 having the argument g. No constant term 
is present in (p/a 0 ) cos^t (3/2) e 0 . Thus we must have 
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After the constants c x and c 2 are determined, we obtain the perturbations of the mean anomaly by 
integrating Equation 69 or 73. The perturbations of the radius vector are obtained from the formula 
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which is convenient for the use of the iteration process. This formula was used by the author in the 
artificial satellites theory (Reference 4). 

Determination of the integration constants associated with the k parameters is done in exactly 
the same way as in Reference 4. We put 
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where[\J (i = 1,2, 3, 4) are the series obtained by the formal integration procedure as applied to 
Equations 37 to 40 and where A and B are the constants of integration. We omit here the development 
leading to the equations for the determination of A and B, because this development can be found in 
Reference 4, and give only the final results: 
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DETERMINATION OF THE POSITION AND VELOCITY VECTORS 
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We determine the position of the satellite with respect to the moving system of coordinates, 
whose x axis is in the orbit plane of the disturbing body and is directed toward the mean position of 
the pericenter; the y axis is in the orbit plane of the disturbing body, normal to the x axis; and the 
z axis is normal, in the positive direction, to the x and y axes. Introducing the matrices 
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For the absolute velocity, but decomposed along the axes of the moving system ( x ,y lZ ), we have 


V - h A 3 (i/r* -n 0 ' -n 0 y't) • A, (I) • A 3 (*- 0 ) 


“ sin f 
cos f + e 
0 


The last equation can be written in the form 


v - hA 3 ( -&>') • A • A, 


A 3 (* -7T o “ n o yt) 


“ sin f 
cos f t e 
0 


Taking the equation 


f + X “ 77 0 - n 0 yt - f 


into account, Equation 83 can be reduced to the form 


V = hAjl-oi'l • A • A 3 Id 


- sin 7 - e sin (x “ w 0 - n 0 yt) 
+ cos T + e cos (x-" 0 -n 0 yt) 
0 


Taking Equations 12 and 13 into account, we have finally: 




1 - e 2 


A 3 { ~ OJ* ) • A • A. ( CU ) 


h_ 

"o 


h; sin f j 


1 “ e 0 2 


A , T \ 1 “ e o 2 ^ 

h„ ( cos f + e o) + 2 T 


(83) 


(84) 


(85) 


0 
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Equation 85 is analogous to Equation 22 of the article published by A. Bailie and R. Bryant (Refer- 
ence 8) on the determination of the velocity of artificial satellites. The transformation from a mov- 
ing system (x.y.z) to the inertial system does not represent any difficulty and therefore is omitted. 

COLLECTION OF FORMULAS 

Input information: 


g 0 . u ' 0 » I 0 . a o- e 0 ' n 0 


/ / • / o 1 


g 0 , a , c , n , i 


Basic arguments: 

g = g 0 + n o t • e' “ %o + n o t 
u) - + n 0 ( y + a - q ) t , 

~ < + n 0 (a + 7 } + y M t , 

6 = (tt 0 ' +n Q y't -a') - a-' . 

We make use of the following standard formulas of the elliptic motion: 


cos (p 


2 e o + 2 


TW 

L 


Jp'(P e o) 


cos py 


sin op 


2 l ^ 1 e o 2 J P (P e o) 

e o L-> p 


sin py 


1 + T e o 2 “ 2e 




j; (p e ») 


cos py 


3 , . V -1 J p ( pe o) 


1 + 2 e o 2 ~ 4 


E 

p=! 


2 cos py 


P 


= 1 + 2 


+ CD 

ZZ j p ( pe °) 


cos py 


p s i 


= 1 + 2 


cos pg 


( pe< ) 

P=1 

If preferred, Cayley’s tables can be used instead of formulas given above. 


i 
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' 1 + l' ) * “ ’ jT COS ( (p + f ' + u + c.P ) , 


s 2 --- (1 + n) 


/ ■•■ a 

“ n _ * -p sin ( cp + f + a' + ) 


: 1 + 1- ) 


A- a 

a n * T 1 


COS ( '• p - f ' t f( : - t/: ' ) 


II + *> * * ^7 


— / * sin ( tp - f ' t (a) - co ' ) , 

+ ( K l " V ) S 1 - 2 \ ^2 S 2 + ( K 4 " s 3 " 2x 3 S 


( 1 + V ) 

a o 


1 d S 

2 <9\, 


+ /V 1 S j ^2 S 2 


_1_ d S 
2 


“S S 1 ’ S 2 


1 Ss 

2 <9\„ 


^3 S 3 " ^4 S 4 




Jl — 

2 Jx~ 


+ X„ s 3 - \ 3 S 4 


2a o n 0 
1 - e 2 


a 2 
a o 


r + 


1 “ e„ 2 - 


a o7 \V 


i! [i -Y 

0 a o 


2a o n o p 


sin <p 


0 yr^7 


i - 


p ' i + 


h 2 

V 


2a o n o 


V 


l - e 0 2 


2 cos <p + 3 e 0 j - xYt Y sin <P 


h 2 

h o 2 


- 1 


sin (p 

1^7 


2a o n o 




P . ^ \ r r ^ d p 

~ cos cp + 2 e Q j + yi-e 0 2 Sine 


h 2 

+ Yr 1 


sin cp d p 
~~ cos 


5 o yrp 


^ o a o P 


1 — p 2 o 2 

1 e 0 a 0 
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n o a 0 P e o sln(p 


1 - e o 2 a ° yro 


3 2 1 2 
2 s 2 P 


m a a 


■ (4 s3 -4 p 2 ) 


rn a"* a ( 35 , 15 2 2 3 • 

= * — — 3 — q * n * + *o~ n ^ 


r* t s “ " 4 " s p + yp 


n = n, + n 2 + o 3 + 


an 


an, 

<?n 2 

<?o 3 

dy 

+ TT 

+ ay 


(90 
^ dp 


20j + 30 2 + 40 3 + 


a r ' dO 
r ’ a ' ‘ dS 


m a a 

+ — • — r * 3s 

o ' 3 r 


0 a 3 a ' / 15 2 3 \ 

’ ~r ’ w s " 2 p j 

n ' a 4 a ' / 35 , 15 2 \ 

VT S -^F sp J 


r = N. • (A = A. 2,3) 


V 1 T 

2_ j n j r-^r(" 0 8«) n (n = 0, 1 , 2, 3, 




dT 


dt = +n 0 y'}' + F, 


dY 

Ht 


H 0 y T + 2 ho • 1 


2 I +F 2 


d ^0 


dt h F 3 ' 


22 


C - -LA 

G i ” 2 h„ 




1 - e 0 2 


_1 Jh_ 

2 h n 


Vi - e o 2 


a 

P a 


1 A 

2 h n 


n - e 0 2 


/a_ # r ' 

1 «Jf)\ r 

V P * 7" 

1 + v (9 S / 


1 h 

a o n o 

/_L 

7' 

1 an \ 


G 4 = 2 h 0 • 

i 1 - e o 2 

\P ' 

a' 

1 + v dS ) 

L + ( 


H i - ^ ( n ° 8z )' 


dk 


J • IT ■ aA Is") ‘ [ + ( S 2 + x 4 2 ) a 2 

- (\, \ 4 + \j* 3 ) <?3 

<%] 

r - • AS) [- K 2 + ^ 3 2 K 

- (\ a \ 4 \ j ) a 3 

+ (\,\ 4 + \ 2 * 3 ) a 4 ] 


( A. j \ 4 + \ 2 \ 3 ) a 3 

( A . 2 A . 4 - A.j \ 3 ) cr 2 J 

! 2+ ^ 2 )-3 

+ (\ 2 \ 3 ) a, 

-(\j\ 4 +\ 2 \ 3 ) a,] ; 

(i = 0,1, 2, 3, 4) , 

an' 


a n g. 

i 

dy n 


Let 


1 1 , V , 

= + n 0 a \ 2 t - 2 ( + \ sin & + cos ‘ sin i « 

1 d .^ ' 

n 0 a\ + H 2 + -2-(+\ 4 cos^“/V 3 sin^) • sin i ' fa ’ 

d \ 3 ^ 

fa = + n o ^4 + ^3 + ~2 ( + cos ^ + A s in 0) * sin i J 

d ^4 1 

r - n 0 tj \ 3 + H 4 + ^ _ (+\ 1 sin ^-\ 2 cos #) * sin i ' 

, [T], W, [A.J (i = 1, 2,3,4) 


d^ 

~dt 


dH 1 

dt 


dir 

dt 


be the series obtained by the formal integration process. 
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1 + c i + hr 


t = c 2 + IT] 


ip = [«p] 


o 3 


[H] = - 3-fj- -ye 0 It] + 2 (a 2 -A 3 


- 3A - -y e.T + 2 (a 2 - A 3 + 


P 

— C°S<J) 


P . - 

— sin qj 


( n o Sz , 

d n p 

- " ! 

dg n a o 

© 

C 

n 

d n p 

n! 

dg n a o 

K Sz ) 

d n p 2 

n! 

dg n a 2 


V “ [ = ] + [T] coscp+^-e 0 ] + [V] sin cp 


B - n Q [T] f— cos f “ ~~ cos cp 


[^3 ( ~PT sin f - — — sin qp 
V d fl “n 


n o y / F 2 p 2 \ n 0 ^ 2 (l+W) 

n r^} U f ~VJ + i - * 2 


dn 0 8 z 


~ n n - 3c 


P - 3 

n o c 2 (-57 cos<p + y e Q 


n ° y p 2 r - 1 

= — T + W J + B 


!\ - e 0 2 a o 


3 o L w oJ 


"o »o 


”0 y p* 

n 2 a 3 + 3 " A i + A 2 cos g 

fl - e 0 2 a o 


— cos (p + y e 0 - /j cos g 
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( 11 ) : 
( 12 ) 


A 1 A 2 e 0 
1 3n ft + 2^n 


W ^ = + T (-£- cos f + 4 e oJ + * a" 0 *“> f 


y (A- W) - y (A + W) v ; 


const, in J f \ 


M)’ * (M " M) '} 
.) - M)’ * (M * [‘J) 


T 1 o + sin 

T I oJ 

+ (11) 

= 0 

l 

1 \ 



2 X o “ sln 

T I o) 

+ (12) 

= 0 


1 1 r , 

\j - sin yl ( + y(A + B) + [.VJ 


[‘J 

M 


\ - cos y I 0 +y(A-B) + [\ 4 J 


+ 1 0 0 
0 + cos cl “ sin a 

Lo + sin a + cos a 


A 3 (a ' - 


+ cos a - s in a 0 

+ sin a + cos a 0 

L 0 0+1 


= A, 


) •A*A 3 (ccj)(1+v) 


r cos f 
7 sin f 

l 0 — I 
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V 


a n n ^ 

- A. I -«'J • A • A, («) 

/T ^7 


- sin f “ 2 ( 1 e o 2 ) ^ 

+ "R^ ( cosT + e o) + T( 1_e o 2 ) T 


*13 
*2 j 
*22 

*23 
*3 1 


= " *, 2 + */ " *3 2 + */ 

= - 2 (Xj ^4 + *1 *2) 

= + 2 (a^j A.j - A.j 

= + 2 (A,A 4 -A | \ a ) ■ 

= “ A,’ + * 2 2 "*3 2 + *4 2 > 

= 2 A.^ + A- 2 A. 3 ) » 

= + 2 (\, ^3 + ^2*4) ’ 

= + 2 (A,A 4 -A 2 * 3 ) • 

= - k* - * 2 2 + *3 2 + *4 2 


CONCLUSION 

A systematic investigation of the orbital stability of satellites under different initial conditions 
can be based on the development given herein. In light of recent results concerning the stability of 
lunar satellites by G. Chebotarev and S. Kirpichnikov (Reference 9), which were obtained by the method 
of numerical integration, it would be of interest to investigate the problem using the development into 
periodic series. We can expect that one of the most influential long period terms in the coordinates 
will be the term with the period equal to the sidereal revolution of the earth around the moon. The 
described method will be programmed by the Theoretical Division, Goddard Space Flight Center, for 
computation on the IBM 7090. We have to expect that the number of terms in the development either 
for planetary satellites (with the exception of the moon) or for lunar satellites will not be large, and 
thus the computation can be done with the limited accuracy. 
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Appendix A 


NOTATIONS 

a 0 the mean semi- major axis of the satellite’s orbit 

a,e,n the osculating elements of the satellite 

e 0 the mean eccentricity of the satellite’s orbit 

f the osculating true anomaly 

f ' the true anomaly of the disturbing body 

e: = g 0 + n Q t the undisturbed mean anomaly of the satellite 

r 0 satellite’s mean anomaly at the epoch 
an 

h 

i 0 the mean value of the mutual inclination of two orbit planes 
-2K the periodic part in 0 - 0' 

l - <* 0 + n 0 t + n 0 6z the perturbed mean anomaly of the satellite 
m' the mass of the disturbing body 
-2N the periodic part in 0 + 0 ' 

n 0 the mean anomalistic mean motion of the satellite (a 0 is defined as n 0 " 2/3 ; the 
gravitational constant is put equal to 1) 

n Q bz the perturbations in the satellite’s mean anomaly 

-2n 0 at the secular part in 0 + 0' 

+2n 0 7;t the secular part in 0 - 0' 

r the position vector of the satellite, r = |r| 

r' the position vector of the disturbing body, |r' | = r' =7' 
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r , f the radius vector and the true anomaly depending upon i , 


S = 

V 

i 

V 

7 

A 

^l» ^2* ^3* ^4 


V 

tt 0 +n 0 yt 

77 o +n oy ,t: 
P. nP 


p, <P 


r cos f - a Q ( cos E - e Q ) 

7 sin? - a Q |/l - e Q 2 sin E 
E - e Q s in E = l 
cos ( r , r ' ) 

the true orbital longitude of the satellite, reckoned from the departure point in 
the satellite's orbit plane 

the true orbital longitude of the disturbing body, reckoned from the departure 
point in the body's orbit plane 

the auxiliary mean anomaly; y is replaced by / at the later stage of computation 

the distance between the satellite and the disturbing body 

the parameters defining the periodic perturbations of the orbit plane, 

1 l 

/Vj - sin'JcosN.Xj - sin2 s inN, 

1 1 

k 3 - cos~2 sin K, k^ - cos 2 cos K 

the perturbations in the radius vector r, 

r = ( 1 + v) 7 

the mean orbital longitude of the pericenter of the satellite 

the mean orbital longitude of the pericenter of the disturbing body 

the radius vector and the true anomaly depending upon the auxiliary mean 
anomaly y, 

p cos p = a Q { cos € ~ e Q ) 
p s in qp = a 0 ^ 1 - e 0 2 sine 
e - e Q sin e = y 

the radius vector and the true anomaly depending upon g, 

p cos qp - a Q (cos7-e 0 ) 
psinqp - a Q |^1 - e Q 2 sine 

e" - e 0 s in 7 - g 
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or' the true orbital longitude of the node of the disturbing body’s orbit plane; the 
fixed ecliptic and equinox serve as a reference 

x the true orbital longitude of the osculating pericenter of the satellite 

the true orbital longitude of the common node, as reckoned from the departure 
point in the orbit plane of the satellite 

0' the true orbital longitude of the common node, as reckoned from the departure 
point in the orbit plane of the disturbing body 

xp 0 , \p' Q the constant parts in 0 and i// 

0 the disturbing function 

the radial component of the disturbing force 

r?Q 

jz the component of the disturbing force normal to the orbit plane 
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NASA- Langley, 1963 G-351 


